#'###########################################
#'                                          #
#'  ------------------------------------    #
#'  |     COIN     __                  |    #
#'  |            >(' )                 |    #
#'  |              )/                  |    #
#'  |             /  `----/            |    #
#'  |             \  ~=- /             |    #
#'  |   ~^~^~^~^~^~^~^~^~^~^~^~^~^~^   |    #
#'  ------------------------------------    #
#'                                          #
#'###########################################



#'#################################################################################################
#'#                                                                                              ##
#'#                                 Fonctions statistiques                                       ##
#'#                                                                                              ##
#'#################################################################################################


#'#################################################################################################
#'#                                                                                              ##
#'#                                 Bootstrap                                                    ##
#'#                                                                                              ##
#'#################################################################################################


#'#
#'# Intervalle de confiance par bootstrap
#'# 
boot.mean.CI <- function (data, type) {
  r <- boot(data, statistic=mean.boot, R=1000);
  ci <- boot.ci(r, type=type)$bca
  return(list(mean = mean(r$t[,1]), lb = ci[4], ub = ci[5]))
}


#'#
#'# Fonction qui dispatche les methodes de calcul de dispersion par bootstrap
#'# 
disp.boot <- function (data, FUN, statBy, key, method="CI", typeCI="bca", ...) {
  tmp <- FUN(data, ...)
  if (!inherits(tmp, "DStat"))
    stop("FUN should return object of class DStat")
  if (method == "CI") {
    res <- tmp[, boot.mean.CI(stat, typeCI), by=statBy]
  } else if (method == "sd") {
    res <- tmp[, {
          m <- mean(stat);
          std <- sd(stat);
          list(mean=m, lb=(m - std), ub=(m + std), sd=std)
        }, by=statBy]
  }
  setkeyv(res, key)
  return (res)
}


#'#
#'# Calcul de la moyenne pour fonction de bootsrap
#'# 
mean.boot <- function(data, indice) {
	mean(data[indice])
}


#'#
#'# Calcul des proportions d'une variable selon un groupe specifie
#'# 
stat.prop <- function(data, countBy, propBy) {
  res <- data[, list(nlocs=nrow(.SD)), by=countBy]
  res <- res[, stat := (nlocs / sum(nlocs))* 100 , by=propBy]
  setattr(res, "class", c("DStat", "data.table", "data.frame"))
  return(res)
}


#'#
#'# Calcul d'un quantile pour les groupes de donnees specifies
#'# 
stat.quantile <- function(data, var, q, quantileBy) {
  if(q > 0) {
    q <- q /100
  }
  expr <- substitute(stat := quantile(var, q))
  res <- data[, eval(expr), by=quantileBy]
  setattr(res, "class", c("DStat", "data.table", "data.frame"))
  return(res)
}




getQuantiles <- function(quantiles, data, quantileBy=quote(list(id, LC)),
    statBy=quote(list(LC)), key=c("LC"), method="sd", type) {
    
    res <- vector("list", length(quantiles))
    for (i in seq_along(quantiles)) {
        tmp <- disp.boot(data, stat.quantile,
            var=distance,
            quantileBy=quantileBy,
            q=quantiles[i],
            statBy=statBy,
            key=key, method=method)
        tmp[, quantile := quantiles[i]]
        tmp[, type := type]
        res[[i]] <- tmp
    }
    return (rbindlist(res))
}

#'#
#'# Creation d'une table de proportions par classe de localisation, session
#'# et visibilite pour tests de compositions statiques
#'# 
prop.table.byLC <- function (data) {
	props <- stat.prop(data, countBy=quote(list(visibility, LC, session, id)),
			propBy=quote(list(visibility, session, id)))
	props[, uid := paste(visibility, session, id, sep="_")]
	setkey(props, uid, LC)
	props <- props[CJ(unique(uid), unique(LC))]
	props <- props[, as.list(setattr(stat, 'names', levels(LC))), by=list(uid)]
	
	ids <- unlist(strsplit(props[, uid], "_"))
	props[, visibility := ids[seq(1,length(ids), by = 3)]]
  props[, session := ids[seq(2,length(ids), by = 3)]]
	props[, id := ids[seq(3,length(ids), by = 3)]]
	
	props$uid <- NULL
	
	return (props)
}



counts.table.byLC <- function (data) {
    counts <- data[, list(nlocs=nrow(.SD)), by=list(visibility, LC, session, id)]
    counts[, uid := paste(visibility, session, id, sep="_")]
    setkey(counts, uid, LC)
    counts <- counts[CJ(unique(uid), unique(LC))]
    counts <- counts[, as.list(setattr(nlocs, 'names', levels(LC))), by=list(uid)]
    
    ids <- unlist(strsplit(counts[, uid], "_"))
    counts[, visibility := ids[seq(1,length(ids), by = 3)]]
    counts[, session := ids[seq(2,length(ids), by = 3)]]
    counts[, id := ids[seq(3,length(ids), by = 3)]]
    
    counts$uid <- NULL
    
    return (counts)
}



#'#
#'# Methode pour standardiser les donnees
#'# 
standardize <- function(data) {
  return((data-mean(data))/sd(data))
}

#'#
#'# Comparaison 2 a 2 de lignes dans une table par le test du Chi carre
#'# 
compare.chisq <- function (table, threshold=0.05) {
  n <- dim(table)[1]
  nres <- factorial(n) / (2 * factorial(n - 2))
  result <- data.frame(var1=character(nres), var2=character(nres), 
      x.squared=numeric(nres), df=numeric(nres), 
      p.value=numeric(nres), stringsAsFactors=FALSE)
  cpt <- 1
  for (i in 1:n) {
    if (i < n) {
      for (j in (i + 1):n){
        tmp <- table[c(i, j), ]
        res <- chisq.test(tmp)
        result[cpt, c("var1", "var2")] <- c(rownames(table)[i], rownames(table)[j])
        result[cpt, c("x.squared", "df", "p.value")] <-  c(round(res$statistic, 2), 
            res$parameter, 
            signif(res$p.value, 2))
        cpt <- cpt + 1
      }
    }
  }
  return (result)
}


#'#
#'# Sequence logarithmique
#'# 
logseq <- function(from, to, length.out) {
  exp(seq(log(from), log(to), length.out = length.out))
} 


#'#
#'# Retourne la proportion de localisation situees dans des intervalles de distance
#'# 
prop.between <- function(FUN, props){
	res <- numeric (length(props) + 1)
	previous <- 0
	for (i in 1:length(props)) {
		tmp <- FUN(props[i])
		res[i] <- tmp - previous
		previous <- tmp
	}
	res [i + 1] <- 1 - sum(res)
	return (res * 100)
}


#'#
#'# Renvoie une table avec les proportions de localisations situees dans les
#'# intervalles fournis par Argos
#'# 
locs.error.prop <- function (data, locsclasses=c(3, 2, 1, 0, "A", "B")) {
	cumuldist <- list()
	for (i in locsclasses) {
		cumuldist[i] <- list(ecdf(data[i][, distance]))
	}
	
	localisations.props <- do.call(rbind, lapply(cumuldist, prop.between, props=c(250, 500, 1500)))
	colnames(localisations.props) <- c("0-250", "250-500", "500-1500", ">1500")
	
	return (localisations.props)
}


#'#
#'# Renvoie les frequences de localisation des tests statiques
#'# 
static.freqs <- function (data) {
	return (disp.boot(data, FUN = stat.prop,
					countBy=quote(list(visibility, LC, id)),
					propBy=quote(list(visibility, id)),
					statBy=quote(list(visibility, LC)),
					key=c("visibility","LC"))) 
}


#'#
#'# Renvoie les rayons d'erreur des tests statiques
#'# 
static.radius <- function(data, quantile=68) {
	return(disp.boot(data, stat.quantile, 
					statBy=quote(list(visibility, LC)),
					key=c("visibility","LC"),
					var=distance,
					quantileBy=quote(list(visibility, id, LC)),
					q=quantile))
}


#'#
#'# Renvoie les frequences de localisation des tests mobiles
#'# 
mobile.freqs <- function (data) {
	return (disp.boot(data, FUN = stat.prop,
					countBy=quote(list(LC, id)),
					propBy=quote(list(id)),
					statBy=quote(list(LC)),
					key=c("LC"))) 
}


#'#
#'# Renvoie les rayons d'erreur des tests mobiles
#'# 
mobile.radius <- function(data, quantile=68) {
	return(disp.boot(data, stat.quantile, 
					statBy=quote(list(LC)),
					key=c("LC"),
					var=distance,
					quantileBy=quote(list(id, LC)),
					q=quantile))
}


#'#
#'# Compare l'efficacite des filtres pour les tests statiques
#'# 
compare.filters.desc.static <- function (x, y) {
	compare.filters.desc(x, y, countBy=quote(list(visibility, LC)),
			propBy=quote(list(visibility)), keys=c("visibility", "LC"))
}


#'#
#'# Compare l'efficacite des filtres pour les tests mobiles
#'# 
compare.filters.desc.mobile <- function (x, y) {
	compare.filters.desc(x, y, countBy=quote(list(LC)),
			propBy=NULL, keys=c("LC"))
}


#'#
#'# Compare le resultat de 2 filtres : nbre de locs supprimees, proportion de locs,
#'# quantite de localisations gagnees ou perdues
#'# 
compare.filters.desc <- function (x, y, countBy, propBy=NULL, keys) {
	props.x <- stat.prop(x, countBy=countBy,
			propBy=propBy)
	props.y <- stat.prop(y, countBy=countBy,
			propBy=propBy)
	
	setkeyv(props.x, keys)
	setkeyv(props.y, keys)
	
	all <- merge(props.x, props.y, by=keys, all.x=TRUE)
	
	replaceall.na.dt(all)
	
	all[, gain := stat.y - stat.x]
	all[, pertelocs := nlocs.x - nlocs.y]
	all[, proplocs := round((pertelocs/nlocs.x) * 100, 2)]
	
	return (all)
}


#'#
#'# Calcule l'effect size pour une anova en fournissant une data.table
#'# la variable et le facteur consideres
#'# 
effect.size.anova <- function(data, variable, category) {
	mexpr <- substitute(mean(variable, na.rm=TRUE))
	nexpr <- substitute(length(variable))
	vexpr <- substitute(var(variable, na.rm=TRUE))
	nomexpr <- substitute(list(nom=(length(variable)/N) * ((mean(variable, na.rm=TRUE) - u)^2)))
	sel <- substitute (!is.na(variable))
	
	u <- data[, eval(mexpr)]
	N <- data[eval(sel), eval(nexpr)]
	v <- data[, eval(vexpr)]
	noms <- data[, eval(nomexpr), by=category]
	f <- noms[, sqrt(sum(nom, na.rm=TRUE)/v)]
	
	return(f)
}


#'#
#'# Moyenne et IC95% pour echnatillon aleatoire du deuxieme degre
#'# 
meanICl2 <- function(dat, var, byG2) {
    
    #' Taille Groupe 1 (deploiements)
    n <- length(levels(dat[[byG2]]))
    #' Taille Groupe 2 (iterations)
    m <- nrow(dat)/n
    #' Ntot Groupe2
    nm <- n * m
    
    #' Moyenne Groupe2
    expr <- parse(text=paste("list(mean=mean(",var,"))", sep=""))
    yi <- dat[, eval(expr), by=byG2]
    #' Moyenne generale
    y <- mean(yi$mean)
    
    #' Somme carres Groupe 2
    expr2 <- parse(text=paste("list(ss=sum((",var," - mean(",var,"))^2))", sep="")) 
    ss2 <- dat[, eval(expr2), by=byG2]
    #' Variance Groupe 2
    s2 <- sum(ss2$ss) / (n*(m - 1))
    #' Variance Groupe 1
    s1 <- sum((yi$mean - y)^2)/(n - 1)
    
    #' Variance Generale
    vg <- s1/n + s2/nm
    #' Erreur type
    se <- sqrt(vg)
    #' Intervalle de confiance (95%)
    ic <- 1.96 * se
    
    return (list(mean=y, var=vg, se=se, li=y-ic, ls=y+ic))
}


#'#
#'# Moyenne et IC95% pour echnatillon aleatoire du troisieme degre
#'# 
meanICl3 <- function(dat, var, byG3, byG2) {
    if (!is.factor(dat[[byG2]])) {
        dat[[byG2]] <- as.factor(dat[[byG2]])
    }
    if (!is.factor(dat[[byG3]])) {
        dat[[byG3]] <- as.factor(dat[[byG3]])
    }
    #' Taille Groupe 1 (deploiements)
    n <- length(levels(dat[[byG2]]))
    #' Taille Groupe 2 (iterations simulation)
    m <- length(levels(dat[[byG3]]))
    #' Ntot Groupe 2
    nm <- n * m
    #' Taille Groupe 3 (iterations HR)
    k <- nrow(dat)/nm
    #' Ntot Groupe 3
    nmk <- n*m*k
    
    #' Moyenne Groupe 3
    expr <- parse(text=paste("list(mean=mean(",var,"))", sep=""))
    yij <- dat[, eval(expr), by=c(byG3, byG2)]
    #' Moyenne Groupe 2
    yi <- dat[, eval(expr), by=byG2]
    #' Moyenne Generale
    y <- mean(yi$mean)
    
    #' Variance Groupe 1
    s1 <- sum((yi$mean - y)^2)/(n - 1)
    #' Variance Groupe 2
    s2 <- sum((yij$mean - yi$mean)^2)/(n*(m - 1))
    #' Somme carres Groupe 3
    expr2 <- parse(text=paste("list(ss=sum((",var," - mean(",var,"))^2))", sep="")) 
    ss3 <- dat[, eval(expr2), by=c(byG3, byG2)]
    #' Variance groupe 3
    s3 <- sum(ss3$ss) / (n*m*(k - 1))
    
    #' Variance Generale
    vg <- s1/n + s2/nm + s3/nmk
    #' Erreur type
    se <- sqrt(vg)
    #' Intervalle de confiance (95%)
    ic <- 1.96*se
    
    return (list(mean=y, var=vg, se=se, li=y-ic, ls=y+ic))
    
}